% This script computes tidal power spectra for all the geomagnetic
% stations using two approaches  1) pmtm 2) robustfit



S = dir('/data/backup/mnair/obs_mag_data/update_2011/secular_removed/*.nc');
a = (datenum(2011,1,1)- datenum(1995,1,1))*24*60;
t = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2010,12,31,23,59,30));
periods = [4 4.8 6 8 11.967236 12  12.421 12.6583 23.934472 24 25.891];

x = [cos(2*pi*t/( periods(1)/24))'       sin(2*pi*t/( periods(1)/24))'];

for j = 2:length(periods),
    
    x = [x cos(2*pi*t/( periods(j)/24))'       sin(2*pi*t/( periods(j)/24))'];
    
end;
%for nstn = 1:141,
for nstn = 141:-1:1,
    
t = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2010,12,31,23,59,30));

    
ncid = netcdf.open(['/data/backup/mnair/obs_mag_data/update_2011/secular_removed/' ... 
    S(nstn).name],'NOWRITE');


X_ID = netcdf.inqVarID(ncid,'Magnetic_Field_X');
Y_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Y');
Z_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Z');
z_data = netcdf.getVar(ncid,Z_ID,0,a);

latitude  = double(netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'geospatial_lat'));
longitude = double(netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'geospatial_lon'));

z_data = double(z_data)/10;
z_data(z_data == 99999.9 ) = NaN;
L = isnan(z_data);
k = find(L==0);

z_data(L) = 0;

netcdf.close(ncid);

if k(end) - k(1) > 4000000,
    st = k(1);
    en = k(1) + 4000000 - 1;
else
    st = k(1);
    en = k(end);
end;

try
    
% [P, F] = pmtm(z_data(st:en), 4, (en-st) + 1, 1/60);
% 
% LL = 1./(F*3600) > 10 & 1./(F*3600) < 30;
%  
% P = P(LL);
% F = F(LL);
 
display(S(nstn).name);
 
catch me    
end;


%robust fit
%all data
clear k LL

if sum(~L) > 4000000,
    L(st+4000000:end) = 1;
end;

% %robust fit
% %all data
% 
% spectra_all = robustfit(x(~L,:), z_data(~L));
% 
% % calculate the amplitudes
% 
% spectra_data_all_rob = abs(complex(spectra_all(2:2:end),spectra_all(3:2:end)));


% %least square solution
% 
% spectra_all_ls = x(~L,:)\z_data(~L);
% 
% spectra_data_all_ls = abs(complex(spectra_all_ls(1:2:end),spectra_all_ls(2:2:end)));
% 

% only day-time data
t = t + (longitude/15)/24; %shift the axis to the local time
%Get the LT

LL = L | (t-floor(t)  > 18/24 & t-floor(t)  < 6/24)';

%robust fit of the night data
spectra_day_rob = robustfit(x(~LL,:), z_data(~LL));
spectra_data_day_rob = abs(complex(spectra_day_rob(2:2:end),spectra_day_rob(3:2:end)));

% only nightime data

%Get the LT

LL = L | (t-floor(t)  < 18/24 & t-floor(t)  > 6/24)';

%robust fit of the night data
spectra_night_rob = robustfit(x(~LL,:), z_data(~LL));
spectra_data_night_rob = abs(complex(spectra_night_rob(2:2:end),spectra_night_rob(3:2:end)));




%ls fit of the night data
% spectra_night_ls = x(~L,:)\z_data(~L);
% spectra_data_night_ls = abs(complex(spectra_night_ls(1:2:end),spectra_night_ls(2:2:end)));

% % spectra using only M2 component
% x = [cos(2*pi*t/( periods(7)/24))'       sin(2*pi*t/( periods(7)/24))'];
% %ls fit of the night data for M2 component
% spectra_night_ls_m2 = x(~L,:)\z_data(~L);
% spectra_data_night_ls_m2 = abs(complex(spectra_night_ls_m2(1:2:end),spectra_night_ls_m2(2:2:end)));

% sweeping in frequency for M2 component LS
% display('Sweep Started');
% prd = 12.415;
% clear spectra_data_ls_m2_array
% for i = 1:10,
%     
% x_m2 = [ones([1,length(t)])' cos(2*pi*t/( prd/24))'  sin(2*pi*t/(prd/24))'    ... %M2
%     ];  
% spectra_night_ls_m2 = x_m2(~L,:)\z_data(~L);
% spectra_data_ls_m2_array(i) = abs(complex(spectra_night_ls_m2(2),spectra_night_ls_m2(3)));
% prd = prd + 0.001;
% end;
% display('Sweep Ended');


% sweeping in frequency for M2 component ROB TOO SLOW
prd = 12.415;
clear spectra_data_rob_m2_n

for i = 1:100,
    
x_m2 = [ cos(2*pi*t/( prd/24))'  sin(2*pi*t/(prd/24))'    ... %M2
    ];  

spectra_night_rob_m2 = robustfit(x_m2(~LL,:),z_data(~LL));
spectra_data_rob_m2_n(i) = abs(complex(spectra_night_rob_m2(2),spectra_night_rob_m2(3)));
prd = prd + 0.0001;
end;

eval(['save /home/mnair/projects/tides/' S(nstn).name(1:end-3) '_ROB_DAY_NIGHT_7_years_SWEEP' ...
    ' spectra_day_rob spectra_data_day_rob spectra_night_rob spectra_data_night_rob spectra_data_rob_m2_n st en latitude longitude periods']);

clear z_data L LL spectra_day_rob spectra_data_day_rob spectra_data_ls_m2_array spectra_data_night_ls spectra_night_ls spectra_data_all_ls spectra_all_ls spectra_data_rob_m2_n;

end;

